library(spdep)
library(MCMCpack)
library(mvtnorm)


CAR.Model <- function(parm, X, y, W, beta.mean=rep(0,times=dim(X)[2]),beta.var=diag(dim(X)[2])*100,sigma.shape=.001,sigma.rate=.001,delta.shape=.1,delta.rate=.1,ridge=.97)
     {
     ### Parameters
     N<-dim(X)[1]
     beta <- parm[1:dim(X)[2]]
     sigma <- exp(parm[dim(X)[2]+1])
     delta <- exp(parm[dim(X)[2]+2])
     phi.vec<- parm[(dim(X)[2]+3):(dim(X)[2]+2+N)]
     Q<-length(phi.vec)

	 ### Log of Prior Densities
     sigma.prior <- log(dinvgamma(sigma, shape=sigma.shape, scale=1/sigma.rate))
     delta.prior <- log(dinvgamma(delta, shape=delta.shape, scale=1/delta.rate))
     log.mu.prior <- sum(log(dmvnorm(beta, mean=beta.mean, sigma = beta.var)))

	### Log-Likelihood
     mu <- tcrossprod(X, t(beta))
     B<-diag(apply(W,1,sum))-ridge*W
     LL <- -(N/2)*log(sigma)-(Q/2)*log(delta)-((t(y-mu-phi.vec)%*%(y-mu-phi.vec))/(2*sigma))-((t(phi.vec)%*%B%*%phi.vec)/(2*delta))      
	
     ### Log-Posterior
     LP <- LL + sigma.prior + delta.prior + log.mu.prior
     return(LP)
     }
     
CAR.Proper<-function(X, y, W, mcmc=1000, burnin=round(.1*mcmc), beta.mean=rep(0,times=dim(X)[2]),beta.var=diag(dim(X)[2])*100,sigma.shape=.001,sigma.rate=.001,delta.shape=.1,delta.rate=.1,ridge=.97, initials=c(rep(.1,dim(X)[2]),1,1,rep(.1,dim(X)[1]))){
	model<-MCMCmetrop1R(CAR.Model, theta.init=initials, y=y, X=X, W=W, beta.mean=beta.mean,beta.var=beta.var, thin=1, mcmc=mcmc, burnin=burnin, verbose=0,logfun=TRUE)
	invisible(model)
	}